module grid
   !
   !  SUBROUTINES
   !
   !  1) set_grid
   !  2) get_uniform_grid
   !  3) get_backward_GP_grid
   !  4) get_GP_ratio
   !
   implicit none

contains

   !****************************************************************************

   subroutine set_grid(folder_output, sim_id, nmbr, kg, nx, ny, a1, rt, x, y) ! Last 2 are output
      implicit none
      character(3), intent(in) :: nmbr
      character(200), intent(in) :: folder_output
      character(200), intent(in) :: sim_id
      integer, intent(in) :: kg ! Kind of grid (1 = uniform)
      integer, intent(in) :: nx
      integer, intent(in) :: ny
      real(8), intent(in) :: a1
      real(8), intent(in) :: rt
      real(8), intent(out) :: x(nx*ny)
      real(8), intent(out) :: y(nx*ny)
      !
      integer :: i, j, np

      ! Read from a file the boundary of the grid
      open(10, file = trim(folder_output) // trim(sim_id) // "-boundary" // &
         nmbr // ".dat")

      ! Reading lower boundary
      j = 1
      read(10,*)
      do i = 1, nx-1
         np = nx*(j-1) + i
         read(10,*) x(np), y(np)
      end do

      ! Reading upper boundary
      read(10,*)
      j = ny - 1
      do i = 1, nx-1
         np = nx*(j-1) + i
         read(10,*) x(np), y(np)
      end do

      if (kg == 1) then

         call get_uniform_grid(nx, ny, x, y)

      else if (kg == 2) then

         call get_backward_GP_grid(a1, nx, ny, x, y)

      else if (kg == 3) then

         call get_power_function_grid(nx, ny, a1, x, y)

      else if (kg == 4) then

         call get_backward_GP_modified_grid(a1, rt, nx, ny, x, y)

      else if (kg == 5) then

         call get_backward_GP_improved_grid(a1, rt, nx, ny, x, y)

      else

         write(*,*) 'set_grid: Type of grid unknown'
         stop

      end if

   end subroutine set_grid

   !============================================================================

   subroutine get_power_function_grid(nx, ny, a1, x, y)
      implicit none
      integer, intent(in) :: nx
      integer, intent(in) :: ny
      real(8), intent(in) :: a1
      real(8), intent(inout) :: x(nx*ny)
      real(8), intent(inout) :: y(nx*ny)

      integer :: i, j, k, np, nps, npi, npf, iaux
      real(8) :: xi, xf, yi, yf, dx, dy
      real(8) :: cte1, cte2, cg, yaux(ny-1), xaux(ny-1)

      iaux = int(a1)
      cte1 = a1 - dble(iaux)
      cte2 = dble(iaux)

      do i = 1, nx-1
         j = 1
         npi = nx*(j-1) + i
         xi = x(npi)
         yi = y(npi)

         j = ny-1
         npf = nx*(j-1) + i
         xf = x(npf)
         yf = y(npf)

         dx = (xf - xi)/(ny - 2)
         dy = (yf - yi)/(ny - 2)

         do j = 2, ny-2
            np = nx*(j-1) + i
            nps = np - nx
            x(np) = x(nps) + dx
            y(np) = y(nps) + dy
         end do

         do j = 1, ny-1
            xaux(j) = (j - 1)*dy/y(npf)
            yaux(j) = 1.0d0/(cte1 + xaux(j)**cte2)
         end do

         do k = 1, 100
            cg = y(npf)/(0.5*yaux(1) + sum(yaux(2:ny-2)) + 0.5*yaux(ny-1))
            do j = 2, ny-2
               np = npi + (j - 1)*nx
               xaux(j) = xaux(j-1) + dy/y(npf)
               yaux(j) = 1.0d0/(cte1 + xaux(j)**cte2)
               dy = 0.5*cg*(yaux(j-1) + yaux(j))
               y(np) = y(np-nx) + dy
            end do
         end do
      end do

   end subroutine get_power_function_grid

   !============================================================================

   subroutine get_uniform_grid(nx, ny, x, y) ! Last 2 are inoutput
      implicit none
      integer, intent(in) :: nx
      integer, intent(in) :: ny
      real(8), intent(inout) :: x(nx*ny)
      real(8), intent(inout) :: y(nx*ny)
      !
      integer :: i, j, np, nps
      real(8) :: xi, xf, yi, yf, dx, dy

      do i = 1, nx-1

         j = 1
         np = nx*(j-1) + i
         xi = x(np)
         yi = y(np)

         j = ny-1
         np = nx*(j-1) + i
         xf = x(np)
         yf = y(np)

         dx = (xf-xi)/(ny-2)
         dy = (yf-yi)/(ny-2)

         do j = 2, ny-2
            np = nx*(j-1) + i
            nps = np - nx
            x(np) = x(nps) + dx
            y(np) = y(nps) + dy
         end do
      end do

   end subroutine get_uniform_grid

   !****************************************************************************

   subroutine get_backward_GP_grid(a1, nx, ny, x, y) ! Last 2 are inoutput
      implicit none
      integer, intent(in) :: nx
      integer, intent(in) :: ny
      real(8), intent(in) :: a1
      real(8), intent(inout) :: x(nx*ny)
      real(8), intent(inout) :: y(nx*ny)
      !
      integer :: i, j, np, npn
      real(8) :: q, xi, xf, yi, yf

      do i = 1, nx-1

         j = 1
         np = nx*(j-1) + i
         xi = x(np)
         yi = y(np)

         j = ny-1
         np = nx*(j-1) + i
         xf = x(np)
         yf = y(np)

         call get_GP_ratio(ny-2, (yf-yi)/a1, q)

         do j = ny-2, 2, -1

            np = nx*(j-1) + i
            npn = np + nx

            x(np) = x(npn)

            y(np) = y(npn) - a1*q**(ny-j-2)

         end do
      end do

   end subroutine get_backward_GP_grid

   !============================================================================

   subroutine get_backward_GP_modified_grid(a1, rt, nx, ny, x, y) ! Last 2 are inoutput
      implicit none
      integer, intent(in) :: nx
      integer, intent(in) :: ny
      real(8), intent(in) :: a1
      real(8), intent(in) :: rt
      real(8), intent(inout) :: x(nx*ny)
      real(8), intent(inout) :: y(nx*ny)

      integer :: i, j, np, npn
      real(8) :: q, xi, xf, yi, yf, aux

      do i = 1, nx-1

         j = 1
         np = nx*(j - 1) + i
         xi = x(np)
         yi = y(np)

         j = ny - 1
         np = nx*(j - 1) + i
         xf = x(np)
         yf = y(np)

         aux = (yf/rt)*a1

         call get_GP_ratio(ny-2, (yf-yi)/aux, q)

         do j = ny-2, 2, -1

            np = nx*(j-1) + i
            npn = np + nx

            x(np) = x(npn)

            y(np) = y(npn) - aux*q**(ny-j-2)

         end do

      end do

   end subroutine get_backward_GP_modified_grid

   !============================================================================

   subroutine get_backward_GP_improved_grid(a1, rt, nx, ny, x, y) ! Last 2 are inoutput
      implicit none
      integer, intent(in) :: nx
      integer, intent(in) :: ny
      real(8), intent(in) :: a1
      real(8), intent(in) :: rt
      real(8), intent(inout) :: x(nx*ny)
      real(8), intent(inout) :: y(nx*ny)

      integer :: i, j, np, npn
      real(8) :: q, xi, xf, yi, yf, aux

      do i = 1, nx-1

         j = 1
         np = nx*(j - 1) + i
         xi = x(np)
         yi = y(np)

         j = ny - 1
         np = nx*(j - 1) + i
         xf = x(np)
         yf = y(np)

         aux = (yf/rt)*a1

         call get_GP_ratio((ny-2)/2, (yf-yi)/aux, q)

         do j = ny-2, 2, -1

            np = nx*(j-1) + i
            npn = np + nx

            x(np) = x(npn)

            if (aux*q**(ny-j-2) < y(npn)/j) then

               y(np) = y(npn) - aux*q**(ny-j-2)

            else

               y(np) = y(npn) - y(npn)/j

            end if

         end do

      end do

   end subroutine get_backward_GP_improved_grid

   !============================================================================

   subroutine get_gp_ratio(n, r, q)
      implicit none
      integer, intent(in) :: n   !< number of partitions
      real(8), intent(in) :: r   !< l/a1
      real(8), intent(out) :: q  !< q

      ! Parameters

      integer :: nit = 1000   ! Maximum number of iteractions
      real(8) :: tol = 1.d-15 ! Tolerance

      ! Inner variables

      integer :: i   ! Dummy index
      real(8) :: qi  ! inital value of q
      real(8) :: qf  ! final value of q
      real(8) :: qm  ! mean value of q

      if (r < n) then

         qi = 0.1d0

         qf = 1.d0 - 1.d-15

      else

         qi = 1.d0 + 1.d-15

         qf = 10.d0

      end if

      do i = 1, nit

         qm = 0.5d0*qi + 0.5d0*qf

         if (0.d0 < f(qi)*f(qm)) then

            qi = qm

         else

            qf = qm

         end if

         if (abs(qf-qi) < tol) exit

      end do


      if (i == nit) then

         write(*,*) "get_gp_ratio: Maximum number of iteractions was exceeded."

         stop

      end if

      q = qm

   contains

      real(8) function f(q)
         implicit none
         real(8), intent(in) :: q

         f = q**n + r*(1.d0 - q) - 1.d0

      end function f

   end subroutine get_gp_ratio

   !============================================================================

end module grid
